R Notebook

library(fpca)

Loading required package: sm
package ‘sm’ was built under R version 3.4.4Package 'sm', version 2.2-5.6: type help(sm) for summary information
Loading required package: splines

# options
# display as decimal, not scientific notation
options(scipen = 999)

# data for 2013 - 2016
county_poverty <- read.csv(file = 'county_poverty_hist.csv', header = TRUE, sep = ',')
county_poverty$geo_sumlevel <- gsub('[[:punct:]]\\s', '', str_extract(county_poverty$geo_name, '\\,\\s[A-Z]{2}'))
county_poverty$geo_name <- gsub('\\,\\s[A-Z]{2}', '', county_poverty$geo_name)
county_poverty <- county_poverty[-3]

# naming scheme:
# pYEAR = number of people in poverty, for the given year
# prYEAR = poverty rate, for the given year
# popYEAR = total population, for the given year
names(county_poverty) <- c('county', 'state', 
                           'p2013', 'p2014', 'p2015', 
                           'pr2013', 'pr2014', 'pr2015', 
                           'pr2016', 'p2016', 'pop2013', 
                           'pop2014', 'pop2015', 'pop2016')

county_poverty <- county_poverty %>% 
  na.omit() %>%
  select(sort(current_vars()))

# remove Puerto Rico
county_poverty <- county_poverty %>% 
  filter(state != 'PR')

county_poverty <- county_poverty %>% 
  select(sort(current_vars()))

# economic/unemployment data
# https://apps.bea.gov/regional/
ec <- read.csv('CAINC4__ALL_STATES_1969_2017.csv', header = TRUE, sep = ',')
ec <- ec[, c('GeoFIPS', 'GeoName', 'Description', 'X2013', 'X2014', 'X2015', 'X2016')]
ec$Description <- gsub('[[:punct:]]', '', ec$Description) %>% 
  trimws()
# keep 'Per capita personal income (dollars) 4/', 'Total employment', 'Farm proprietors' income', 'Nonfarm proprietors' income'
# clean data
ec <- ec %>% 
  filter(ec$Description == 'Per capita personal income dollars 4' | 
                      ec$Description == 'Total employment' | 
                      ec$Description == 'Farm proprietors income' | 
                      ec$Description == 'Nonfarm proprietors income')
ec$Description[ec$Description == 'Per capita personal income dollars 4'] <- 'income'
ec$Description[ec$Description == 'Total employment'] <- 'empl_count'
ec$Description[ec$Description == 'Farm proprietors income'] <- 'farm_income'
ec$Description[ec$Description == 'Nonfarm proprietors income'] <- 'nonfarm_income'
ec <- dcast(melt(ec, 
                 id.vars = c('GeoFIPS', 'GeoName', 'Description'), 
                 measure.vars = c('X2013', 'X2014', 'X2015', 'X2016')), 
            GeoFIPS + GeoName + variable ~ Description) %>% 
  mutate(variable = gsub('X', '', variable))

attributes are not identical across measure variables; they will be dropped

ec[, 1:2] <- lapply(ec[, 1:2], as.character)
ec[, c(1, 3:7)] <- lapply(ec[, c(1, 3:7)], as.numeric)

NAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercion

ec$p_farm_income <- ec$farm_income / ec$nonfarm_income
ec[is.na(ec)] <- 0
# combine economic data with land area
area <- read_excel('LND01.xls')[,c(2,4)]

Error in read_fun(path = path, sheet = sheet, limits = limits, shim = shim,  : 
  path[1]="LND01.xls": No such file or directory

# create tall data set
cp1 <- gather(county_poverty[c(1:5, 14)], year, poverty_count, p2013:p2016, factor_key = TRUE) %>% 
  mutate(year = str_extract(year, '\\d{4}'))
cp2 <- gather(county_poverty[c(1, 6:9, 14)], year, total_count, pop2013:pop2016, factor_key = TRUE) %>% 
  mutate(year = str_extract(year, '\\d{4}'))
cp <- merge(cp1, cp2, by = c('county', 'state', 'year')) %>% 
  mutate(poverty_rate = poverty_count / total_count,
         year = year %>% as.numeric())
cp$county <- gsub('\\sCounty', '', cp$county)
cp$county <- gsub('\\sParish', '', cp$county)
cp$county <- gsub('\\scity', '\\sCity', cp$county)
cp$county <- gsub('[[:punct:]]', '', cp$county)
cp$county <- cp$county %>% 
  trimws()
# county data
county_ec <- ec[!(ec$GeoFIPS %% 1000 == 0), ]
county_ec <- subset(county_ec, select = -c(GeoFIPS))
county_ec$state <- gsub('[[:punct:]]\\s', '', str_extract(county_ec$GeoName, '\\,\\s[A-Z]{2}'))
county_ec$GeoName <- gsub('\\,\\s[A-Z]{2}', '', county_ec$GeoName)
names(county_ec) <- c('county', 'year', 'empl_count', 'farm_income', 'income', 'nonfarm_income', 'p_farm_income', 'area', 'state')
county_ec$county <- gsub('[[:punct:]]', '', county_ec$county)
county_ec$county <- gsub('\\sIndependent\\sCity', '\\sCity', county_ec$county)
county_ec$county <- gsub('\\sStaunton\\sWaynesboro', '', county_ec$county)

# manual fix for counties with area = 0
# scrape areas from Wikipedia, as Census data was inadequate
zeroarea <- county_ec[which(county_ec$area == 0), -c(2:7)] %>% 
  unique()
zeroarea$wiki_link <- gsub(' City and Borough', '', zeroarea$county)
zeroarea$wiki_link <- gsub(' Municipality', '', zeroarea$wiki_link)
zeroarea$wiki_link <- gsub('[A-Z]{2}', '', 
                           paste('https://en.wikipedia.org/wiki/', 
                                 gsub('\\s', '_', zeroarea$wiki_link), ',_', 
                                 gsub('AK', 'Alaska', zeroarea$state), 
                                 gsub('CO', 'Colorado', zeroarea$state), 
                                 sep = ''))
for(i in 1:nrow(zeroarea)) {
  zeroarea$area[i] <- gsub(',', '', read_html(zeroarea$wiki_link[i]) %>% 
                        html_nodes(xpath = '//th[contains(., "Area")]//parent::tr/following-sibling::tr[1]/td') %>% 
                        html_text() %>% 
                        str_extract('[\\d\\.,]+'))[1] %>% 
  as.numeric()
  county_ec$area[which(county_ec$county == zeroarea$county[i])] <- zeroarea$area[which(zeroarea$county == zeroarea$county[i])]
}

county_data <- merge(cp, county_ec, by = c('county', 'state', 'year')) %>% 
  mutate(empl_rate = empl_count / total_count, 
         pop_dens = total_count / area)
# 'county_data' to wide
st <- list()
for (i in 2013:2016) {
  df <- county_data %>% 
    filter(year == i)
  df <- df[, -3]
  df_names <- names(df)
  df_names[3] <- paste('p', i %>% as.character(), sep = '')
  df_names[4] <- paste('pop', i %>% as.character(), sep = '')
  df_names[5] <- paste('pr', i %>% as.character(), sep = '')
  df_names[6] <- paste('ec', i %>% as.character(), sep = '')
  df_names[7] <- paste('f_inc', i %>% as.character(), sep = '')
  df_names[8] <- paste('inc', i %>% as.character(), sep = '')
  df_names[9] <- paste('nf_inc', i %>% as.character(), sep = '')
  df_names[10] <- paste('pf_inc', i %>% as.character(), sep = '')
  df_names[12] <- paste('er', i %>% as.character(), sep = '')
  df_names[13] <- paste('pop_dens', i %>% as.character(), sep = '')
  names(df) <- df_names
  st[[i - 2012]] <- df
}
county_data_wide <- merge(st[[1]], st[[2]], by = c('state', 'county', 'area'))
county_data_wide <- merge(county_data_wide, st[[3]], by = c('state', 'county', 'area'))
county_data_wide <- merge(county_data_wide, st[[4]], by = c('state', 'county', 'area'))

# data set of poverty, by state
state_poverty <- table(county_data_wide$state) %>% 
  data.frame()
names(state_poverty) <- c('state', 'n_counties')
# 51 = number of states, including Washington D.C.
# create a data table of summary stats per state
sp <- list()
for(i in 1:51) {
  state_stats <- data.frame(summary(county_data_wide[which(county_data_wide$state == state_poverty$state[i]), ]))[13:258, 2:3] %>% 
  mutate(stat = gsub('\\s*\\:', '', str_extract(Freq, '.*\\:')), 
         Freq = gsub(':', '', str_extract(Freq, '\\:.+')) %>% as.numeric()) %>% 
  reshape(idvar = 'Var2', timevar = 'stat', direction = 'wide') %>%
  mutate(Year = str_extract(Var2, '\\d{4}'), 
         Var2 = str_extract(Var2, '[a-z_]{1,7}'), 
         state = state_poverty$state[i])
  sp[[i]] <- state_stats
}
# clean data, rename variables and reorganize column order
state_poverty <- merge(state_poverty, do.call(rbind, sp), by = 'state') %>% 
  select(sort(current_vars()))
state_poverty <- state_poverty[, c(8, 10, 7, 9, 1:6)]
state_poverty$Var2[state_poverty$Var2 == 'p'] <- 'poverty_count'
state_poverty$Var2[state_poverty$Var2 == 'pr'] <- 'poverty_rate'
state_poverty$Var2[state_poverty$Var2 == 'pop'] <- 'total_count'
state_poverty$Var2[state_poverty$Var2 == 'ec'] <- 'empl_count'
state_poverty$Var2[state_poverty$Var2 == 'f_inc'] <- 'farm_income'
state_poverty$Var2[state_poverty$Var2 == 'inc'] <- 'income'
state_poverty$Var2[state_poverty$Var2 == 'f_inc'] <- 'farm_income'
state_poverty$Var2[state_poverty$Var2 == 'nf_inc'] <- 'nonfarm_income'
state_poverty$Var2[state_poverty$Var2 == 'pf_inc'] <- 'p_farm_income'
state_poverty$Var2[state_poverty$Var2 == 'er'] <- 'empl_rate'
state_poverty$Var2[state_poverty$Var2 == 'area'] <- 'area'
state_poverty$Var2[state_poverty$Var2 == 'pop_den'] <- 'pop_dens'
names(state_poverty) <- c('state', 'year', 'n_counties', 'countystat', '1st_quantile', '3rd_quantile', 'max', 'mean', 'median', 'minimum')
# replace NAs for 'year'
for(i in 1:nrow(state_poverty)) {
  if(state_poverty$year[i] %>% is.na) {
    state_poverty$year[i] <- state_poverty$year[i + 1]
  }
}
# create variables measuring the standard deviation and variance of county stats, per year
state_variance <- c()
state_sd <- c()
for(i in 1:nrow(state_poverty)) {
  state_variance[i] <- county_data %>% 
    filter(state == state_poverty$state[i] & year == state_poverty$year[i]) %>% 
    select(state_poverty$countystat[i]) %>% 
    var()
  state_sd[i] <- sapply(county_data %>% filter(state == state_poverty$state[i] & year == state_poverty$year[i]) %>% select(state_poverty$countystat[i]), sd)[[1]][1]
}
state_poverty$var <- state_variance
state_poverty$sd <- state_sd
state_poverty[is.na(state_poverty)] <- 0
# merge state stats with statistics on total state population ('state_ec'), total state poverty count, and overall state poverty rate (aggregating 'cp')
# change 'state_ec' names from full words to 2-letter abbreviations
for(i in 1:nrow(state_ec)) {
  if(state_ec$state[i] == 'District of Columbia') {
    state_ec$state[i] <- 'DC'
  } else {
    state_ec$state[i] <- state.abb[grep(state_ec$state[i], state.name)]
  }
}

number of items to replace is not a multiple of replacement lengthnumber of items to replace is not a multiple of replacement lengthnumber of items to replace is not a multiple of replacement lengthnumber of items to replace is not a multiple of replacement length

state_ec <- merge(state_ec, 
                  merge(aggregate(poverty_count ~ state + year, cp, sum), 
                        aggregate(total_count ~ state + year, cp, sum), 
                        by = c('state', 'year')), 
                  by = c('state', 'year')) %>% 
  mutate(poverty_rate = poverty_count / total_count, 
         empl_rate = empl_count / total_count)
state_ec <- merge(state_poverty, state_ec, by = c('state', 'year'))
# rename state-wide variables for clarity
names(state_ec) <- c('state', 'year', 'n_counties', 'countystat', '1st_quantile', '3rd_quantile', 'max', 'mean', 'median', 'minimum', 'var', 'sd', 'state_empl_count', 'state_farm_income', 'state_avg_income', 'state_nonfarm_income', 'state_p_farm_income', 'state_area', 'state_poverty_count', 'state_total_count', 'state_poverty_rate', 'state_empl_rate')

row.names(state_data) <- state_data$state

non-unique values when setting 'row.names': ‘AK’, ‘AL’, ‘AR’, ‘AZ’, ‘CA’, ‘CO’, ‘CT’, ‘DC’, ‘DE’, ‘FL’, ‘GA’, ‘HI’, ‘IA’, ‘ID’, ‘IL’, ‘IN’, ‘KS’, ‘KY’, ‘LA’, ‘MA’, ‘MD’, ‘ME’, ‘MI’, ‘MN’, ‘MO’, ‘MS’, ‘MT’, ‘NC’, ‘ND’, ‘NE’, ‘NH’, ‘NJ’, ‘NM’, ‘NV’, ‘NY’, ‘OH’, ‘OK’, ‘OR’, ‘PA’, ‘RI’, ‘SC’, ‘SD’, ‘TN’, ‘TX’, ‘UT’, ‘VA’, ‘VT’, ‘WA’, ‘WI’, ‘WV’, ‘WY’Error in `row.names<-.data.frame`(`*tmp*`, value = c(1L, 1L, 1L, 1L, 1L,  : 
  duplicate 'row.names' are not allowed

lm(state_vcrime_rate ~ state + year + var, 
   data = state_data %>% filter(countystat == 'poverty_rate')) %>% 
  msummary()

               Estimate  Std. Error t value             Pr(>|t|)    
(Intercept) -0.16600898  0.02987540  -5.557 0.000000121172815565 ***
stateAL     -0.00256400  0.00016705 -15.349 < 0.0000000000000002 ***
stateAR     -0.00185904  0.00020580  -9.033 0.000000000000000738 ***
stateAZ     -0.00282847  0.00017330 -16.321 < 0.0000000000000002 ***
stateCA     -0.00280151  0.00019896 -14.081 < 0.0000000000000002 ***
stateCO     -0.00389012  0.00017447 -22.297 < 0.0000000000000002 ***
stateCT     -0.00434224  0.00032718 -13.272 < 0.0000000000000002 ***
stateDC      0.00709309  0.00037155  19.090 < 0.0000000000000002 ***
stateDE     -0.00146553  0.00036157  -4.053 0.000080609750724398 ***
stateFL     -0.00233992  0.00020156 -11.609 < 0.0000000000000002 ***
stateGA     -0.00348890  0.00018308 -19.056 < 0.0000000000000002 ***
stateHI     -0.00444839  0.00020108 -22.123 < 0.0000000000000002 ***
stateIA     -0.00396951  0.00029232 -13.579 < 0.0000000000000002 ***
stateID     -0.00478225  0.00022306 -21.440 < 0.0000000000000002 ***
stateIL     -0.00296354  0.00022463 -13.193 < 0.0000000000000002 ***
stateIN     -0.00304924  0.00025604 -11.909 < 0.0000000000000002 ***
stateKS     -0.00313147  0.00024636 -12.711 < 0.0000000000000002 ***
stateKY     -0.00520906  0.00018352 -28.384 < 0.0000000000000002 ***
stateLA     -0.00191695  0.00017918 -10.699 < 0.0000000000000002 ***
stateMA     -0.00281371  0.00026630 -10.566 < 0.0000000000000002 ***
stateMD     -0.00233531  0.00018345 -12.730 < 0.0000000000000002 ***
stateME     -0.00553603  0.00030059 -18.417 < 0.0000000000000002 ***
stateMI     -0.00248043  0.00024486 -10.130 < 0.0000000000000002 ***
stateMN     -0.00444102  0.00028297 -15.694 < 0.0000000000000002 ***
stateMO     -0.00219310  0.00020779 -10.555 < 0.0000000000000002 ***
stateMS     -0.00474653  0.00021443 -22.135 < 0.0000000000000002 ***
stateMT     -0.00375591  0.00017411 -21.572 < 0.0000000000000002 ***
stateNC     -0.00346685  0.00020981 -16.524 < 0.0000000000000002 ***
stateND     -0.00459135  0.00017237 -26.637 < 0.0000000000000002 ***
stateNE     -0.00403924  0.00026724 -15.115 < 0.0000000000000002 ***
stateNH     -0.00463935  0.00033665 -13.781 < 0.0000000000000002 ***
stateNJ     -0.00434025  0.00023178 -18.725 < 0.0000000000000002 ***
stateNM     -0.00078217  0.00016708  -4.681 0.000006292790458193 ***
stateNV     -0.00010532  0.00026927  -0.391              0.69625    
stateNY     -0.00301692  0.00025059 -12.039 < 0.0000000000000002 ***
stateOH     -0.00400565  0.00022659 -17.678 < 0.0000000000000002 ***
stateOK     -0.00250750  0.00022000 -11.398 < 0.0000000000000002 ***
stateOR     -0.00418656  0.00028998 -14.437 < 0.0000000000000002 ***
statePA     -0.00356266  0.00029106 -12.240 < 0.0000000000000002 ***
stateRI     -0.00448894  0.00024850 -18.064 < 0.0000000000000002 ***
stateSC     -0.00183783  0.00019820  -9.273 < 0.0000000000000002 ***
stateSD     -0.00445006  0.00053285  -8.351 0.000000000000040245 ***
stateTN     -0.00062988  0.00023868  -2.639              0.00919 ** 
stateTX     -0.00292348  0.00016753 -17.450 < 0.0000000000000002 ***
stateUT     -0.00470278  0.00019850 -23.692 < 0.0000000000000002 ***
stateVA     -0.00509680  0.00017780 -28.666 < 0.0000000000000002 ***
stateVT     -0.00553576  0.00032929 -16.811 < 0.0000000000000002 ***
stateWA     -0.00399023  0.00023101 -17.273 < 0.0000000000000002 ***
stateWI     -0.00389852  0.00026240 -14.857 < 0.0000000000000002 ***
stateWV     -0.00362839  0.00022946 -15.812 < 0.0000000000000002 ***
stateWY     -0.00481323  0.00021980 -21.898 < 0.0000000000000002 ***
year         0.00008573  0.00001484   5.775 0.000000042337838656 ***
var          0.17040076  0.08546137   1.994              0.04796 *  

Residual standard error: 0.0002361 on 151 degrees of freedom
Multiple R-squared:  0.9892,    Adjusted R-squared:  0.9855 
F-statistic: 267.2 on 52 and 151 DF,  p-value: < 0.00000000000000022

lm(state_vcrime_rate ~ state + year + var, 
   data = state_data %>% filter(countystat == 'pop_dens')) %>% 
  msummary()

                    Estimate       Std. Error t value             Pr(>|t|)    
(Intercept) -0.1827872107075  0.0303324801240  -6.026     0.00000001231859 ***
stateAL     -0.0025493182041  0.0001668002380 -15.284 < 0.0000000000000002 ***
stateAR     -0.0020979024549  0.0001667972137 -12.578 < 0.0000000000000002 ***
stateAZ     -0.0029184502459  0.0001668019413 -17.497 < 0.0000000000000002 ***
stateCA     -0.0028741075426  0.0001807304400 -15.903 < 0.0000000000000002 ***
stateCO     -0.0039084943962  0.0001715695923 -22.781 < 0.0000000000000002 ***
stateCT     -0.0048580721635  0.0001682370770 -28.876 < 0.0000000000000002 ***
stateDC      0.0064310304350  0.0001667963389  38.556 < 0.0000000000000002 ***
stateDE     -0.0020356735023  0.0001701650369 -11.963 < 0.0000000000000002 ***
stateFL     -0.0025321098431  0.0001675660907 -15.111 < 0.0000000000000002 ***
stateGA     -0.0033014425537  0.0001677913113 -19.676 < 0.0000000000000002 ***
stateHI     -0.0046580265936  0.0001669320526 -27.904 < 0.0000000000000002 ***
stateIA     -0.0044460794373  0.0001667988342 -26.655 < 0.0000000000000002 ***
stateID     -0.0050759510004  0.0001667974490 -30.432 < 0.0000000000000002 ***
stateIL     -0.0032158361676  0.0001683788492 -19.099 < 0.0000000000000002 ***
stateIN     -0.0034169488730  0.0001670618208 -20.453 < 0.0000000000000002 ***
stateKS     -0.0034863016927  0.0001668255622 -20.898 < 0.0000000000000002 ***
stateKY     -0.0050451230154  0.0001668970430 -30.229 < 0.0000000000000002 ***
stateLA     -0.0017758906576  0.0001668862211 -10.641 < 0.0000000000000002 ***
stateMA     -0.0026129656517  0.0003419663795  -7.641     0.00000000000231 ***
stateMD     -0.0020062377339  0.0002870198856  -6.990     0.00000000008269 ***
stateME     -0.0060332541499  0.0001667973205 -36.171 < 0.0000000000000002 ***
stateMI     -0.0028024760819  0.0001676667651 -16.715 < 0.0000000000000002 ***
stateMN     -0.0048562286115  0.0001679446942 -28.916 < 0.0000000000000002 ***
stateMO     -0.0023798550820  0.0001693199909 -14.055 < 0.0000000000000002 ***
stateMS     -0.0044775551126  0.0001667965858 -26.844 < 0.0000000000000002 ***
stateMT     -0.0038547370663  0.0001667963284 -23.110 < 0.0000000000000002 ***
stateNC     -0.0037039233242  0.0001669848897 -22.181 < 0.0000000000000002 ***
stateND     -0.0045058887798  0.0001667963239 -27.014 < 0.0000000000000002 ***
stateNE     -0.0044474679715  0.0001668400463 -26.657 < 0.0000000000000002 ***
stateNH     -0.0052165157168  0.0001668195873 -31.270 < 0.0000000000000002 ***
stateNJ     -0.0029459452911  0.0008497418100  -3.467             0.000686 ***
stateNM     -0.0007657653601  0.0001668009087  -4.591     0.00000922818479 ***
stateNV     -0.0005257684625  0.0001667967761  -3.152             0.001955 ** 
stateNY      0.0112834420635  0.0071310049495   1.582             0.115672    
stateOH     -0.0042770322678  0.0001676181947 -25.517 < 0.0000000000000002 ***
stateOK     -0.0027859216144  0.0001668337819 -16.699 < 0.0000000000000002 ***
stateOR     -0.0046375660160  0.0001671314123 -27.748 < 0.0000000000000002 ***
statePA     -0.0035689816481  0.0002824186807 -12.637 < 0.0000000000000002 ***
stateRI     -0.0047882652724  0.0001700080362 -28.165 < 0.0000000000000002 ***
stateSC     -0.0020469628469  0.0001668071027 -12.271 < 0.0000000000000002 ***
stateSD     -0.0034410995171  0.0001667963070 -20.631 < 0.0000000000000002 ***
stateTN     -0.0009597260154  0.0001668712505  -5.751     0.00000004759584 ***
stateTX     -0.0029277776899  0.0001671974513 -17.511 < 0.0000000000000002 ***
stateUT     -0.0048991442243  0.0001670196539 -29.333 < 0.0000000000000002 ***
stateVA     -0.0046485052382  0.0003234171580 -14.373 < 0.0000000000000002 ***
stateVT     -0.0061011334634  0.0001667965308 -36.578 < 0.0000000000000002 ***
stateWA     -0.0042988167742  0.0001668642171 -25.762 < 0.0000000000000002 ***
stateWI     -0.0042972375386  0.0001668136639 -25.761 < 0.0000000000000002 ***
stateWV     -0.0039400637857  0.0001667998460 -23.622 < 0.0000000000000002 ***
stateWY     -0.0050985128235  0.0001667963351 -30.567 < 0.0000000000000002 ***
year         0.0000943869796  0.0000150569724   6.269     0.00000000362896 ***
var         -0.0000000002637  0.0000000001281  -2.058             0.041290 *  

Residual standard error: 0.0002359 on 151 degrees of freedom
Multiple R-squared:  0.9893,    Adjusted R-squared:  0.9856 
F-statistic: 267.6 on 52 and 151 DF,  p-value: < 0.00000000000000022

lm(state_vcrime_rate ~ state + year + mean, 
   data = state_data %>% filter(countystat == 'p_farm_income')) %>% 
  msummary()

               Estimate  Std. Error t value             Pr(>|t|)    
(Intercept) -0.12103252  0.03431532  -3.527             0.000557 ***
stateAL     -0.00244393  0.00016927 -14.438 < 0.0000000000000002 ***
stateAR     -0.00197214  0.00017094 -11.537 < 0.0000000000000002 ***
stateAZ     -0.00289953  0.00016505 -17.568 < 0.0000000000000002 ***
stateCA     -0.00294148  0.00016706 -17.607 < 0.0000000000000002 ***
stateCO     -0.00393034  0.00016628 -23.636 < 0.0000000000000002 ***
stateCT     -0.00490279  0.00016487 -29.738 < 0.0000000000000002 ***
stateDC      0.00643088  0.00016487  39.007 < 0.0000000000000002 ***
stateDE     -0.00206329  0.00016553 -12.465 < 0.0000000000000002 ***
stateFL     -0.00249554  0.00016672 -14.969 < 0.0000000000000002 ***
stateGA     -0.00324939  0.00016792 -19.351 < 0.0000000000000002 ***
stateHI     -0.00467270  0.00016487 -28.343 < 0.0000000000000002 ***
stateIA     -0.00425524  0.00017857 -23.830 < 0.0000000000000002 ***
stateID     -0.00487300  0.00018018 -27.046 < 0.0000000000000002 ***
stateIL     -0.00314827  0.00016987 -18.534 < 0.0000000000000002 ***
stateIN     -0.00335816  0.00016720 -20.085 < 0.0000000000000002 ***
stateKS     -0.00327491  0.00018218 -17.976 < 0.0000000000000002 ***
stateKY     -0.00500224  0.00016602 -30.131 < 0.0000000000000002 ***
stateLA     -0.00173208  0.00016603 -10.433 < 0.0000000000000002 ***
stateMA     -0.00322817  0.00016487 -19.581 < 0.0000000000000002 ***
stateMD     -0.00245627  0.00016523 -14.866 < 0.0000000000000002 ***
stateME     -0.00602590  0.00016489 -36.544 < 0.0000000000000002 ***
stateMI     -0.00282840  0.00016490 -17.152 < 0.0000000000000002 ***
stateMN     -0.00477015  0.00017090 -27.913 < 0.0000000000000002 ***
stateMO     -0.00236022  0.00016728 -14.109 < 0.0000000000000002 ***
stateMS     -0.00443108  0.00016572 -26.739 < 0.0000000000000002 ***
stateMT     -0.00367647  0.00017663 -20.814 < 0.0000000000000002 ***
stateNC     -0.00364337  0.00016712 -21.801 < 0.0000000000000002 ***
stateND     -0.00436628  0.00017217 -25.360 < 0.0000000000000002 ***
stateNE     -0.00386418  0.00026730 -14.456 < 0.0000000000000002 ***
stateNH     -0.00522206  0.00016487 -31.675 < 0.0000000000000002 ***
stateNJ     -0.00465974  0.00016487 -28.264 < 0.0000000000000002 ***
stateNM     -0.00066467  0.00016894  -3.934             0.000127 ***
stateNV     -0.00055989  0.00016529  -3.387             0.000900 ***
stateNY     -0.00336834  0.00016504 -20.409 < 0.0000000000000002 ***
stateOH     -0.00429453  0.00016497 -26.032 < 0.0000000000000002 ***
stateOK     -0.00266827  0.00017076 -15.626 < 0.0000000000000002 ***
stateOR     -0.00457832  0.00016737 -27.355 < 0.0000000000000002 ***
statePA     -0.00402305  0.00016495 -24.389 < 0.0000000000000002 ***
stateRI     -0.00485617  0.00016487 -29.455 < 0.0000000000000002 ***
stateSC     -0.00203865  0.00016492 -12.361 < 0.0000000000000002 ***
stateSD     -0.00303679  0.00021883 -13.877 < 0.0000000000000002 ***
stateTN     -0.00096544  0.00016487  -5.856         0.0000000286 ***
stateTX     -0.00284826  0.00016892 -16.862 < 0.0000000000000002 ***
stateUT     -0.00480114  0.00016994 -28.253 < 0.0000000000000002 ***
stateVA     -0.00520559  0.00016493 -31.562 < 0.0000000000000002 ***
stateVT     -0.00607852  0.00016507 -36.824 < 0.0000000000000002 ***
stateWA     -0.00424245  0.00016654 -25.474 < 0.0000000000000002 ***
stateWI     -0.00425670  0.00016566 -25.696 < 0.0000000000000002 ***
stateWV     -0.00394624  0.00016487 -23.935 < 0.0000000000000002 ***
stateWY     -0.00507286  0.00016511 -30.723 < 0.0000000000000002 ***
year         0.00006373  0.00001703   3.741             0.000259 ***
mean        -0.00020643  0.00007347  -2.810             0.005616 ** 

Residual standard error: 0.0002332 on 151 degrees of freedom
Multiple R-squared:  0.9895,    Adjusted R-squared:  0.9859 
F-statistic:   274 on 52 and 151 DF,  p-value: < 0.00000000000000022

lm(state_vcrime_rate ~ state + year + var, 
   data = state_data %>% filter(countystat == 'empl_rate')) %>% 
  msummary()

               Estimate  Std. Error t value             Pr(>|t|)    
(Intercept) -0.16584128  0.03046629  -5.443  0.00000020690236819 ***
stateAL     -0.00232649  0.00028211  -8.247  0.00000000000007374 ***
stateAR     -0.00187305  0.00028276  -6.624  0.00000000057725935 ***
stateAZ     -0.00269917  0.00027963  -9.653 < 0.0000000000000002 ***
stateCA     -0.00281504  0.00026403 -10.662 < 0.0000000000000002 ***
stateCO     -0.00383766  0.00022849 -16.796 < 0.0000000000000002 ***
stateCT     -0.00467380  0.00028559 -16.365 < 0.0000000000000002 ***
stateDC      0.00668593  0.00030639  21.822 < 0.0000000000000002 ***
stateDE     -0.00187060  0.00028960  -6.459  0.00000000136363296 ***
stateFL     -0.00235305  0.00027167  -8.662  0.00000000000000661 ***
stateGA     -0.00318462  0.00022905 -13.903 < 0.0000000000000002 ***
stateHI     -0.00442819  0.00029723 -14.898 < 0.0000000000000002 ***
stateIA     -0.00422026  0.00028415 -14.852 < 0.0000000000000002 ***
stateID     -0.00520768  0.00021352 -24.390 < 0.0000000000000002 ***
stateIL     -0.00304324  0.00027795 -10.949 < 0.0000000000000002 ***
stateIN     -0.00322055  0.00027460 -11.728 < 0.0000000000000002 ***
stateKS     -0.00328264  0.00027010 -12.153 < 0.0000000000000002 ***
stateKY     -0.00485795  0.00026157 -18.572 < 0.0000000000000002 ***
stateLA     -0.00160218  0.00025089  -6.386  0.00000000199158404 ***
stateMA     -0.00307112  0.00023035 -13.332 < 0.0000000000000002 ***
stateMD     -0.00226145  0.00028240  -8.008  0.00000000000028970 ***
stateME     -0.00580699  0.00028394 -20.451 < 0.0000000000000002 ***
stateMI     -0.00261719  0.00027827  -9.405 < 0.0000000000000002 ***
stateMN     -0.00467744  0.00027726 -16.870 < 0.0000000000000002 ***
stateMO     -0.00222256  0.00027575  -8.060  0.00000000000021527 ***
stateMS     -0.00427042  0.00026828 -15.918 < 0.0000000000000002 ***
stateMT     -0.00364065  0.00027312 -13.330 < 0.0000000000000002 ***
stateNC     -0.00351146  0.00026909 -13.049 < 0.0000000000000002 ***
stateND     -0.00435704  0.00022522 -19.345 < 0.0000000000000002 ***
stateNE     -0.00422841  0.00028351 -14.914 < 0.0000000000000002 ***
stateNH     -0.00499346  0.00028503 -17.519 < 0.0000000000000002 ***
stateNJ     -0.00443768  0.00028050 -15.820 < 0.0000000000000002 ***
stateNM     -0.00056434  0.00026533  -2.127              0.03505 *  
stateNV     -0.00122579  0.00072234  -1.697              0.09176 .  
stateNY     -0.00323640  0.00022823 -14.180 < 0.0000000000000002 ***
stateOH     -0.00409035  0.00027859 -14.682 < 0.0000000000000002 ***
stateOK     -0.00258715  0.00026694  -9.692 < 0.0000000000000002 ***
stateOR     -0.00443097  0.00028468 -15.565 < 0.0000000000000002 ***
statePA     -0.00382159  0.00027513 -13.890 < 0.0000000000000002 ***
stateRI     -0.00461594  0.00029420 -15.690 < 0.0000000000000002 ***
stateSC     -0.00183411  0.00027537  -6.660  0.00000000047726822 ***
stateSD     -0.00324702  0.00025773 -12.598 < 0.0000000000000002 ***
stateTN     -0.00075813  0.00027152  -2.792              0.00591 ** 
stateTX     -0.00277209  0.00024685 -11.230 < 0.0000000000000002 ***
stateUT     -0.00472267  0.00025786 -18.315 < 0.0000000000000002 ***
stateVA     -0.00503090  0.00025307 -19.879 < 0.0000000000000002 ***
stateVT     -0.00590757  0.00025777 -22.918 < 0.0000000000000002 ***
stateWA     -0.00408172  0.00028349 -14.398 < 0.0000000000000002 ***
stateWI     -0.00407603  0.00028290 -14.408 < 0.0000000000000002 ***
stateWV     -0.00373564  0.00026742 -13.969 < 0.0000000000000002 ***
stateWY     -0.00490828  0.00025467 -19.273 < 0.0000000000000002 ***
year         0.00008585  0.00001514   5.669  0.00000007092442572 ***
var          0.00243152  0.00244261   0.995              0.32111    

Residual standard error: 0.0002384 on 151 degrees of freedom
Multiple R-squared:  0.989, Adjusted R-squared:  0.9853 
F-statistic:   262 on 52 and 151 DF,  p-value: < 0.00000000000000022

state_data %>% 
  filter(countystat == 'poverty_rate') %>% 
  select(var) %>% 
  unique() %>% 
  ggplot(aes(x = var)) + 
  geom_density()

plot(as.party(rpart(state_vcrime_rate ~ ., data = state_data)))

plot(as.party(rpart(state_vcrime_rate ~ var + state_avg_income, data = state_data %>% filter(countystat == 'income'))))

state_data_pca <- state_data_pca_df[, -2] %>% 
  prcomp(scale = TRUE)
-1 * state_data_pca$rotation[, 1:4] %>% round(3)

                       PC1    PC2    PC3    PC4
year                -0.002 -0.049  0.190 -0.593
n_counties          -0.246  0.111 -0.228  0.179
state_total_count   -0.419 -0.037  0.120  0.051
state_area          -0.088 -0.024 -0.051 -0.037
state_pop_dens       0.068 -0.529 -0.006  0.106
state_avg_income     0.037 -0.333  0.509  0.110
state_p_farm_income  0.076  0.141 -0.177  0.673
state_empl_count    -0.417 -0.047  0.138  0.062
state_empl_rate      0.131 -0.468  0.119  0.273
state_poverty_count -0.426 -0.039  0.050  0.034
state_poverty_rate  -0.137 -0.057 -0.550 -0.196
violent_crime       -0.419 -0.074  0.066  0.014
state_vcrime_rate   -0.018 -0.474 -0.245 -0.109
property_crime      -0.421 -0.049  0.022  0.047
state_pcrime_rate   -0.012 -0.338 -0.449 -0.066

summary(state_data_pca)

Importance of components:
                          PC1    PC2    PC3    PC4     PC5    PC6    PC7     PC8     PC9    PC10    PC11    PC12    PC13    PC14    PC15
Standard deviation     2.3236 1.7514 1.4704 1.0927 1.04345 0.9558 0.6842 0.55525 0.44079 0.34091 0.24546 0.13071 0.08364 0.05586 0.02395
Proportion of Variance 0.3599 0.2045 0.1442 0.0796 0.07259 0.0609 0.0312 0.02055 0.01295 0.00775 0.00402 0.00114 0.00047 0.00021 0.00004
Cumulative Proportion  0.3599 0.5644 0.7086 0.7882 0.86077 0.9217 0.9529 0.97343 0.98638 0.99413 0.99815 0.99929 0.99975 0.99996 1.00000

aggregate(pca_data, by = list(pca_data$rowname), FUN = mean) %>% 
  ggplot(aes(x = -PC1, y = -PC2)) + 
  geom_text(aes(label = Group.1), size = 3)

argument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NA

pca_data_alt <- aggregate(pca_data, by = list(pca_data$rowname), FUN = mean)

argument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NAargument is not numeric or logical: returning NA

rownames(pca_data_alt) <- pca_data_alt$Group.1
pca_data_alt <- pca_data_alt[, 3:17]
pca_data_alt %>% 
  dist() %>% 
  hclust() %>% 
  plot(cex = 0.8)

pca_data2 %>% 
  dist() %>% 
  hclust() %>% 
  plot(cex = 0.8, labels = pca_data2$rowname)

NAs introduced by coercion

# visualization
df <- aggregate(state_poverty_rate ~ state, data = state_data, mean)
df <- filter(df, 
             state != "DC") 
mycolors_True <- brewer.pal(9, "GnBu")
states <- tolower(rownames(USArrests))
map_50S <- map_data("state")
ggplot(data = df,
       aes(fill = state_poverty_rate)) +
  labs(x = "Longitude", y = "Latitude", title = "Average Poverty Rate in States") +
  geom_map(aes(map_id = states), map = map_50S) +
  expand_limits(x = map_50S$long, y = map_50S$lat) +
  theme_bw() + 
  scale_fill_gradientn(name = "Poverty Rate", colours = mycolors_True,
                       na.value = "grey") 

aggregate(state_pcrime_rate ~ state, data = state_data, mean) %>% 
  ggplot(aes(x = state, 
             y = state_pcrime_rate, 
             colour = state, 
             group = state)) +
  geom_label(aes(label = state), check_overlap = TRUE)+
  labs(x = "State", 
       y = "Average State Property Crime Rate Per state (including DC)", 
       title = "Average State Property Crime Rate Per state (including DC)")

Ignoring unknown parameters: check_overlap